PrecSezione 4.3: Trasformata di Fourier di sequenze SuCapitolo 4: Campionamento, quantizzazione ed elaborazione numerica Section 4.5: Filtraggio numerico via DFTSegue

## 4.4Trasformata discreta di Fourier↓

Disponendo di una sequenza di lunghezza finita composta da N valori xn, n= 0, 1, ...,N− 1, corrispondenti a campioni di un segnale x(t) prelevati con ritmo fc=(1)/(Tc), si indica come Discrete Fourier Transform (DFT↓) la nuova sequenza [72] [72]La (6.22↓) può essere fatta discendere dalla (6.20↑) vincolando f ad assumere i soli valori discreti f=(m)/(N)(1)/(Tc), e limitando l’indice della sommatoria ad un insieme finito di campioni.
(6.22) Xm=N− 1⎲⎳m= 0xne−j2π(m)/(N)n
univocamente definita per m= 0, 1, ...,N− 1, e che costituisce una approssimazione [73] [73] Una prima fonte di approssimazione deriva dall’operazione di finestratura legata all’uso di un numero finito di campioni, operando quindi su xw(t)=x(t)w(tc) anziché su x(t). Per analizzare le altre fonti di approssimazione, iniziamo a scrivere l’espressione di Xw(f)= ℱ{xw(t)} per f=(m)/(N)fc:
Xw⎛⎝f=(m)/(N)fc⎞⎠ = (N− 1)Tc⌠⌡0x(t)e−j2π(m)/(N)fctdt ≃ N− 1⎲⎳n= 0xn⋅(N− 1)Tc⌠⌡0sinc(fc(t−nTc))e−j2π(n)/(N)fctdt
in cui la seconda eguaglianza utilizza l’interpolazione cardinale x(t)=∑∞n=  − ∞xn⋅sinc(fc(t−nTc)) fornita dalla (6.15↑), ed introduce una seconda fonte di approssimazione legata all’intervallo finito di variazione per m (infatti, benchè l’integrale abbia estensione limitata, i valori di x(t) che cadono entro tale estensione, dovrebbero dipendere da tutti i suoi campioni). L’ultimo integrale è a sua volta una approssimazione (a causa degli estremi di integrazione limitati, e peggiore per i sinc centrati in prossimità dei confini della finestra) della trasformata (calcolata in f=(m)/(N)fc) di sinc(fc(t−nTc)), pari quest’ultima a Tcrectfc(f)e−j2πfnTc, che quando valutata per f=(m)/(N)fc, fornisce il risultato
Xw⎛⎝f=(m)/(N)fc⎞⎠≃TcN− 1⎲⎳n= 0xne−j2π(m)/(N)n
per valori |m|≤(N)/(2), a causa della estensione limitata (in frequenza) di rectfc(f). E’ però facile verificare che Xw⎛⎝(n)/(N)fc⎞⎠ è periodica in m con periodo N, cosicché i valori assunti per m=(N)/(2)+ 1,(N)/(2)+ 2, … sono uguali a quelli per m=  −(N)/(2)+ 1,  −(N)/(2)+ 2, …. del campionamento in frequenza della trasformata X(f)= ℱ{x(t)}, calcolata per f=(m)/(N)fc, e moltiplicata per fc:
(6.23) Xm≃fcX⎛⎝f=(m)/(N)fc⎞⎠
Notiamo subito che la (6.22↑) è valida per qualsiasi m, ed ha un andamento periodico con periodo N, a cui corrisponde la frequenza fc=(1)/(Tc), in accordo con la periodicità in frequenza evidenziata per X•(f) (vedi (6.20↑) e fig. 4.2↑); per questo motivo, qualora il segnale originario x(t) contenga componenti a frequenze maggiori di (fc)/(2), gli Xm con indici prossimi ad (N)/(2) presenteranno errore di aliasing [74] [74]Come osservato al § 4.1.1↑, lo spettro X•(f) di un segnale campionato a frequenza fc è costituito dalle repliche del segnale originario, distanziate di multipli di fc: X•(f)=∑∞n=  − ∞X(f−nfc), e coincide con X(f) per −fc⁄ 2 <f<fc⁄ 2, se X(f) è limitata in banda tra ±W ed fc≥ 2W. Al contrario, se fc< 2W, allora le repliche X(f−nfc) si sovrappongono, e la (6.23↑) si riscrive come Xm≃fcX•⎛⎝f=(m)/(N)fc⎞⎠.. Notiamo inoltre come i valori Xm siano tutti relativi a frequenze ≥ 0, ma nel caso di una sequenza xn reale la proprietà di simmetria coniugata, associata alla periodicità in frequenza, rende il risultato interessante solo per indici m≤(N)/(2), dato che successivamente si trovano valori coniugati a quelli della prima metà. Notiamo infine che la (6.22↑) può essere espressa in forma matriciale: ad esempio, per N= 4 si ottiene
(6.24) ⎡⎢⎢⎢⎢⎢⎣ X0 X1 X2 X3 ⎤⎥⎥⎥⎥⎥⎦=⎡⎢⎢⎢⎢⎢⎣ 1 1 1 1 1 e− j(π)/(2) e−jπ e−j(3π)/(2) 1 e−jπ e−j2π e−j3π 1 e−j(3π)/(2) e−j3π e−j(9π)/(2) ⎤⎥⎥⎥⎥⎥⎦⎡⎢⎢⎢⎢⎢⎣ x0 x1 x2 x3 ⎤⎥⎥⎥⎥⎥⎦
da cui notiamo la proprietà di simmetria per la matrice dei coefficienti.

a - modulo della DFT di una sinusoide a 10 Hz, fc= 100 Hz, con finestre di diversa lunghezza |  |
b - modulo della ℱ-trasformata di una sinusoide a 10 Hz, con finestre di diversa lunghezza |  |

Figura 4.16 Confronto tra DFT ed ℱ-trasformata con uguale estensione temporale
Esempio Allo scopo di concretizzare le differenze tra la trasformata di Fourier ed i valori forniti dalla dft, in fig. 4.16↑-a sono riportati i valori |Xn| per la dft di una sinusoide a 10 Hz, adottando due diverse finestre di analisi (vedi § 3.9.3↑), prelevando alla medesima frequenza di campionamento (100 Hz) un numero variabile di campioni (mostrato in figura), e ponendo i rimanenti a zero, per calcolare in tutti i casi la medesima DFT a 256 punti [75] [75]Il metodo esposto di porre a zero i campioni fino al raggiungimento di una potenza di due, in modo da utilizzare per il calcolo una FFT, è detto zero padding. Il calcolo della DFT su di un numero di punti pari ai campioni di segnale disponibili, non avrebbe dato luogo all’effetto finestra, ma avrebbe fornito in tutti i casi andamenti simili a quello osservabile per 256 punti. Infine, notiamo che nelle figure sono mostrati solo i primi 128 valori, essendo i rimanenti speculari..
Il risultato è quindi confrontato (fig. 4.16↑-b) con quello ottenibile per via analitica calcolando la ℱ-trasformata dello stesso segnale, adottando le medesime finestre temporali, di durata uguale al primo caso. Le curve ottenute nel caso di 80 msec (e 8 campioni!) dipendono da meno di un periodo di segnale, e perciò presentano una componente continua. Aumentando la durata della finestra, l’approssimazione di calcolare una ℱ{} mediante la dft migliora, anche se persiste un ridotto potere di risoluzione spettrale.
Osservazione Probabilmente a questo punto più di un lettore può trovarsi stupito di non avere incontrato due linee spettrali, come sarebbe stato lecito aspettarsi per una sinusoide. La spiegazione è che ciò può accadere, a patto che il numero di campioni N su cui si effettua la dft sia un multiplo intero k del numero di campioni M=T⁄Tc che ricadono entro uno stesso periodo T della sinusoide (con la ovvia condizione che sia M> 2 per rispettare il vincolo fc>2⁄T); in questo caso infatti, l’applicazione della idft (6.25↓) produce una sequenza ancora periodica. Nel caso in esame, essendo T= 100 msec e scegliendo N=  64 e k= 6, il vincolo N=kM=kTfc permette di ottenere fc=(N)/(kT)=(64)/(6⋅10− 1)= 106,6 Hz, ovvero M=Tfc= 10− 1⋅106,6= 10,6 campioni/periodo. La fig. 4.17↓ mostra questo risultato, evidenziando come la riga spettrale si manifesti per m= 6, ossia alla frequenza f=(m)/(N)fc=(6)/(64)106.6 = 10 Hz, mentre la riga presente in m= 58 è in realtà il ripiegamento periodico di quella a −  10 Hz.

Figura 4.17 Sinusoide campionata su di un numero intero di periodi e relativo modulo di dft
Il passaggio dai campioni nel tempo xn a quelli in frequenza Xm è invertibile [76] [76]Sostituendo infatti la (6.22↑) nella (6.25↓), otteniamo
xn=(1)/(N)N− 1⎲⎳m= 0(N− 1⎲⎳k= 0xke−j2π(m)/(N)k)ej2π(m)/(N)n=(1)/(N)N− 1⎲⎳k= 0xkN− 1⎲⎳m= 0ej2π(m)/(N)(k−n)
ma, dato che ∑N− 1m= 0ej2π(m)/(N)(k−n)=⎧⎨⎩ N se k=n 0 altrimenti , allora nella sommatoria esterna sopravvive solo il termine xn, dimostrando l’uguaglianza., ricorrendo alla Inverse Discrete Fourier Trasform (IDFT↓)
(6.25) xn=(1)/(N)N− 1⎲⎳m= 0Xmej2π(m)/(N)n
che per n esterno a [0,N− 1] continua a valere, ed assume valori periodici, coerentemente a quanto accade per lo sviluppo in serie di Fourier. Infatti il legame tra dft e serie di Fourier è molto stretto, in quanto i valori Xm rappresentano una approssimazione [77] [77]La relazione (6.26↓) si dimostra combinando le relazioni (6.7↑) e (6.23↑): Xn≃fcX⎛⎝(n)/(N)fc⎞⎠=fcX⎛⎝(n)/(NTc)⎞⎠=fcX⎛⎝(n)/(T)⎞⎠=fcX(nF)=fcTXSFn=(1)/(Tc)TXSFn=NXSFn dei rispettivi coefficienti della serie di Fourier XSFm=(1)/(T)∫(T)/(2)−(T)/(2)xT(t)e−j2π(m)/(T)tdt, calcolati a partire da un segmento xT(t) estratto da x(t), e moltiplicati per N:
(6.26) Xn≃N⋅XSFn
Per approfondire i risvolti di questo risultato, affrontiamo la sezione successiva.

### 4.4.1Relazione tra DTFT, DFT e trasformata zeta

Così come per i segnali analogici sussiste una relazione (vedi pag. 1↓) tra la trasformata di Fourier e quella di Laplace, così nel contesto delle sequenze, esistono legami tra dtft e trasformata zeta↓, definita quest’ultima come
(6.27) X(z)=∞⎲⎳n=  − ∞xnz−n
che, nel caso in cui converga per |z|= 1, può essere fatta corrispondere alla dtft (6.20↑) semplicemente ponendo z=ejω, ovvero calcolando X(z) sul cerchio unitario↓:
(6.28) X(ejω)=∞⎲⎳n=  − ∞xne−jωn=X(z)|z=ejω=X•(f)|f=(ω)/(2πTc)
Infatti, nelle consuete condizioni in cui gli xn sono i campioni di un segnale x(t) limitato in banda tra ±W e prelevati con ritmo fc≥ 2W, la (6.28↑) effettivamente coincide (per −π≤ω<π) con la X•(f) (eq. (6.20↑), per −(fc)/(2)≤f≤(fc)/(2)) in cui si ponga f=(ω)/(2πTc), mettendo cioè in corrispondenza le frequenze ±(fc)/(2) di X•(f) con le pulsazioni ±π di X(ejω). Al di fuori di tale intervallo, X(ejω) è periodica in ω con periodo 2π, analogamente a ciò che risulta (con periodo fc) per la trasformata di Fourier X•(f) di sequenze; se invece xn è sottocampionata, ossia fc< 2W, anche X(ejω) è affetta da aliasing, così come avviene per X•(f).
Esempio Consideriamo la sequenza xn=⎧⎨⎩ an se  ≥ 0 0altrim. il cui andamento per a= 0.7 è mostrato in fig. 4.18↑, la cui trasformata zeta X(z)=∑∞n=  − ∞anz−n risulta pari a [78] [78]Il risultato si ottiene ricordando che ∑∞n= 0αn=(1)/(1 −α) qualora |α|< 1. X(z)=(1)/(1 −az− 1), ed il cui modulo, dopo aver scritto la variabile complessa z come z=x+jy, è espresso come |X(z)|=1⁄√(⎛⎝(x2−ax+y2)/(x2+y2)⎞⎠2+⎛⎝(ay)/(x2+y2)⎞⎠2).
Figura 4.18 Sequenza xn=an e modulo della relativa trasformata di Fourier a tempo discreto

Facendo ora variare z nell’intervallo [− 2 −j2, 2  +j2], si ottiene per il modulo di X(z) l’andamento mostrato nella figura a lato, in cui è anche raffigurato un cilindro di raggio unitario, la cui intersezione con |X(z)| individua l’andamento di |X(ejω)|=1⁄√1 +a2− 2acosω, ossia la dtft della sequenza an, che a sua volta è mostrata in fig. 4.18↑ per −π<ω<π. |

Se la X(ejω), ottenuta per una sequenza xn aperiodica nel tempo, è campionata in N punti equispaziati e disposti sul cerchio unitario, ossia ponendo ω= 2π(m)/(N) con m= 0, 1, …,N− 1, si ottiene allora una sequenza periodica in frequenza↓ [79] [79]Infatti, e−j2π(m+N)/(N)n=e−j2π(m)/(N)ne−j2πn, ed il secondo termine vale 1 per qualsiasi n. Indichiamo qui ed al prossimo §, una sequenza periodica mediante la tilde .̃.
(6.29) X̃m=∞⎲⎳n=  − ∞xne−j2π(m)/(N)n=X(ejω)|ω=  2π(m)/(N)=X•(f)|f=fc(m)/(N)
che può coincidere con la sequenza Xm ottenuta calcolando la dft (6.22↑) di una sequenza xn, qualora questa abbia una durata limitata ≤N. D’altro canto, è possibile applicare la idft (6.25↑) ad un periodo della sequenza X̃m, ed ottenere quindi una nuova sequenza di valori nel tempo, anch’essa periodica di periodo N, espressa come
(6.30) x̃n=(1)/(N)N− 1⎲⎳m= 0X̃mej2π(m)/(N)n
Infatti, i valori x̃n dipendono da quelli xn=x(t)|t=nTc del segnale originario x(t), campionato agli istanti t=nTc, mediante la relazione [80] [80]Infatti, sostituendo la (6.29↑) in (6.30↑), otteniamo
x̃n=(1)/(N)N− 1⎲⎳k= 0∞⎲⎳h=  − ∞xhe−j2π(k)/(N)hej2π(k)/(N)n
che, scambiando l’ordine delle sommatorie, riscriviamo come
x̃n=∞⎲⎳h=  − ∞xh⎛⎝(1)/(N)N− 1⎲⎳k= 0e−j2π(k)/(N)(h−n)⎞⎠
Dato che (1)/(N)∑N− 1k= 0e−j2π(k)/(N)(h−n)=⎧⎨⎩ 1 se h=n+rN 0 altrimenti , con r intero, si ottiene il risultato (6.31↓).
(6.31) x̃n=∞⎲⎳r=  − ∞xn+rN
e quindi i primi N valori di x̃n coincidono con i campioni di x(t) solo se quest’ultimo ha durata limitata, con estensione minore di NTc, ossia se N è sufficientemente elevato in modo che NTc copra tutta la durata di x(t), e la (6.27↑) si riconduca alla somma di un numero finito di termini. D’altra parte, se x(t) ha durata maggiore di NTc, ovvero X(z) è stata campionata su di un numero di campioni troppo ristretto, allora l’applicazione della IDFT (6.30↑) ad X̃(k) provoca il fenomeno di aliasing temporale↓.
Esempio Nella parte sinistra di fig. 4.19↓viene mostrato il modulo della sequenza Xm ottenuta campionando la X(ejω) dell’esempio precedente, utilizzando 16, 8 o 4 campioni/periodo. Nella parte destra della stessa figura sono quindi rappresentate le corrispondenti sequenze xn ottenute mediante idft. Si può notare che, mentre con 16 campioni/periodo la ricostruzione della sequenza xn=an è piuttosto fedele, con 8 campioni si inizia a verificare il fenomeno di aliasing temporale, che diviene ancor più evidente per 4 campioni/periodo.
Figura 4.19 Aliasing temporale al diminuire della risoluzione del campionamento in frequenza

### 4.4.2Fast Fourier Transform

La sigla fft↓ descrive una classe di algoritmi di calcolo della dft e della sua inversa, caratterizzati dall’uso di un numero molto ridotto di operazioni, rendendo grazie a ciò l’elaborazione numerica dei segnali computatazionalmente praticabile. Analizziamo innanzitutto come il calcolo di ognuno degli N termini Xm della (6.22↑), considerando i valori e−j2π(m)/(N)n precalcolati (vedi (6.24↑)), richieda N moltiplicazioni complesse (equivalenti ognuna ad 4 moltiplicazioni e 2 somme reali) ed N− 1 somme complesse (ognuna pari a 2 somme reali): pertanto una dft richiede N(N(4 + 2)+ 2(N− 1))=N(8N− 2)≃8N2 operazioni. Al contrario, gli algoritmi fft più efficienti riducono il numero di operazioni ad 8Nlog2N: ad esempio, ponendo N= 1024, si ottiene un miglioramento di 23(210)2⁄23210⋅10=210⁄10≃100 volte! Queste prestazioni sono legate all’adozione di un valore di N che sia una potenza di due (ossia N= 2Mcon M intero), ma successivamente sono stati individuati metodi [81] [81]Vedi ad es. `http://it.wikipedia.org/wiki/Trasformata_di_Fourier_veloce` che permettono una efficienza di calcolo comparabile anche per finestre di analisi di lunghezza qualsiasi.

### 4.4.3Relazione tra DFT e DCT

Anche per la dft risulta valida la proprietà di simmetria coniugata (§ 37↑) e quindi, se i valori della sequenza xn di lunghezza N che compare nella (6.22↑) sono reali anziché complessi, allora i coefficienti di dft Xm presentano parte reale pari e parte immaginaria dispari. In particolare, se immaginiamo di estendere la lunghezza della sequenza a 2N punti, ottenuti ribaltando sugli indici negativi la sequenza di partenza come x−n=xn (vedi prima riga di fig. 4.20↓), allora siamo nelle condizioni di sequenza reale pari, che determina una trasformata solo reale (e pari), con parte immaginaria nulla.
Figura 4.20 Estensione pari di sequenza reale
Per arrivare a definire la Discrete Cosine Transform (dct↓) si calcola una dft bilatera sulla sequenza lunga 2N ottenuta traslando quella descritta in modo da renderla effettivamente pari (seconda riga di fig. 4.20↑). Considerando che per segnali reali pari le componenti sinusoidali della base della dft non danno contributi al risultato [82] [82]Scriviamo la (6.22↑) come
Xm = N−1⁄2⎲⎳n’  =  −N+1⁄2xn’  −1⁄2e−j2π(n’)/(2N)m= = N−1⁄2⎲⎳n’  =  −N+1⁄2xn’  −1⁄2cos⎛⎝2π(n’)/(2N)m⎞⎠−jN−1⁄2⎲⎳m’ =  −N+1⁄2xn’  −1⁄2sin⎛⎝2π(n’)/(2N)m⎞⎠= = 2N−1⁄2⎲⎳n’  =1⁄2xn’  −1⁄2cos⎛⎝2π(n’)/(2N)m⎞⎠= 2N− 1⎲⎳n= 0xncos⎛⎝2π(n+1⁄2)/(2N)m⎞⎠= = 2N− 1⎲⎳n= 0xncos⎡⎣(π)/(N)⎛⎝n+(1)/(2)⎞⎠m⎤⎦
in cui xn’ è quella disegnata per seconda in fig. 4.20↑. La quarta eguaglianza tiene conto del fatto che il termine immaginario si annulla, in quanto sommatoria bilatera di una funzione dispari (ottenuta come prodotto di xn’  −1⁄2 pari e sin⎛⎝2π(n’)/(2N)m⎞⎠ dispari), e del fatto che essendo i termini coseno pari, la sommatoria può essere ristretta ai soli indici positivi, raddoppiati. La penultima eguaglianza rappresenta il semplice cambio di variabile n=n’ −1⁄2, mentre l’ultima è (a parte il fattore 2) la definizione della DCT data in (6.32↓)., e adottando un nuovo cambio di variabile, si ottiene in definitiva la formula di calcolo della dct come
(6.32) Xm=N− 1⎲⎳n= 0xncos⎡⎣(π)/(N)⎛⎝n+(1)/(2)⎞⎠m⎤⎦
a cui è associata la trasformazione inversa idct
(6.33) xn=(1)/(2)X0+N− 1⎲⎳m= 1Xmcos⎡⎣(π)/(N)⎛⎝m+(1)/(2)⎞⎠n⎤⎦
La dct verrà usata in questo testo nell’ambito della compressione di immagini (§ 377↓): infatti i valori di luminanza dei pixel in cui si scompone una immagine, sono tutti valori reali.
PrecSezione 4.3: Trasformata di Fourier di sequenze SuCapitolo 4: Campionamento, quantizzazione ed elaborazione numerica Section 4.5: Filtraggio numerico via DFTSegue